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We study theoretically the ac Josephson effect in super- 
conducting planar d-wave junctions. The insulating barrier 
assumed to be present between the two superconductors may 
have arbitrary strength. Many properties of this system de- 
pend on the orientation of the d-wave superconductor: we 
calculate the ac components of the Josephson current. In 
some arrangements there is substantial negative differential 
conductance due to the presence of mid-gap states. We study 
how robust these features are to finite temperature and also 
comment on how the calculated current-voltage curves com- 
pare with experiments. For some other configurations (for 
small barrier strength) we find zero-bias conductance peaks 
due to multiple Andreev reflections through midgap states. 
Moreover, the odd ac components are strongly suppressed and 
even absent in some arrangements. This absence will lead to 
a doubling of the Josephson frequency. All these features are 
due to the d-wave order parameter changing sign when rotated 
90°. Recently, there have been several theoretical reports on 
parallel current in the d-wave case for both the stationary 
Josephson junction and for the normal metal-superconductor 
junction. Also in our case there may appear current density 
parallel to the junction, and we present a few examples when 
this takes place. Finally, we give a fairly complete account 
of the method used and also discuss how numerical calcula- 
tions should be performed in order to produce current-voltage 
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I. INTRODUCTION 

There are now strong reasons to believe that many 
of the high-T c superconductors (HTS) exhibit a super- 
conducting order parameter which has maialy (or com- 
pletely) d a .2_ y 2-wave (d-wave) symmetry.™ The exact 
structure of the order parameter depends on the mate- 
rial. For instance, due to the orthorhombic character of 
YBa 2 Cu307_(5 (YBCO) an s-wave component of the or- 
der parameter is induced B This is not the case for tetrag- 
onal systems like T^F^CuOe+a, where a pure d x i_ y i- 
wave symmetry has been established.^ At any rate, the 
presence of a substantial component of d-wave symmetry 
in the high-T c superconductors cannot be ignored. 

This paper sets out to explore theoretically the impli- 
cations of d-wave symmetry for the ac Josephson effect 
in high-T c superconductors. 



To have the ac Josephson effect, two superconductors 
must be present on both sides of a normal (or insulating) 
region as shown in Fig. [j]. For an anisotropic supercon- 
ductor the order parameter A(9) depends on the angle 
9 of incidence of the quasiparticle approaching the junc- 
tion region. In the case of pure d-wave symmetry we 
have A(9) ~ Aq cos[2(# — a)], where a is the orientation 
of the d-wave order parameter with respect to the insula- 
tor. We will in this paper use the notation d a for d-wave 
superconductor with orientation a. For the isotropic s- 
wave case, A (9) = A s . Another competing candidate 
for the order parameter of HTS has been the anisotropic 
s-wave order parameter: A(9) — A cos 4 [2(9 — a)] + A1 ; 
today, however, it does not seem to be on the short list 
of possible HTS order parameters any longer ac 

The strength of the insulating barrier shown in Fig. |l| 
may vary in our treatment. This makes it possi- 
ble to treat all superconductor-superconductor junctions 
from the ballistic (fully transparent) superconductor- 
normal metal-superconductor (SNS) junction to the 
superconductor-insulator-superconductor (SIS) tunnel 
junction. To describe superconductor-superconductor 
junctions in general we use the notation Si | S2 where Si (2) 
is the superconductor on the left (right) side of the inter- 
face | . The interface is a barrier of arbitrary strength. 

Insulator in short normal region 




Superconductor 1 




Superconductor 2 
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FIG. 1. Layout of the junction. The orientations of the two 
d-wave superconductors are described by the angles ot\ and 
q?2. An electron- like quasiparticle is incident on the junction 
at an angle 9. Short normal regions are introduced between 
the barrier and the superconductors. The strength of the 
barrier can be tuned from the ballistic to the insulating case. 

The ac Josephson effect, from a general point of view, 
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means that for an S|S junction a time independent bias 
voltage V results in a time dependent current. The non- 
stationary current perpendicular to the junction (the IE- 
direction; see Fig. |l|) can be decomposed into frequency 
components 

I x (V,t) = ^2l x , m (V)e im ^ t , (1) 

m 

where m is an integer. The Josephson frequency ujj — 
2eV/h corresponds to the energy gained by a Cooper pair 
when passing through the junction region. Since the re- 
lation I x - m {V) — IxmOO holds, the current is real. 

An order parameter with d^^^-wave symmetry 
changes sign when rotated 90°. This feature has been 
shown to have several implications on transport proper- 
ties of dj;2_ a 2-wave superconductors, some of which have 
been confirmed in experiments: cancellation effects as to 
the dc Josephson currenirTQ and the presence of zero- 
bias conductanee-fjeaks (ZBCP) in various structures of 
HTS j unctions BllB Since in this paper the focus is on the 
ac Josephson effect we will (among other things) discuss 
issues related to ZBCP. 

The ZBCP arises due to the presence of zero-energy 
bound states, so-called midgap states (MGS), at sur- 
faces/interfaces of d x 2_ y 2-w&ve superconductors.EI The 
reason MGS appear is an interplay between normal scat- 
tering (at the surface/interface) and the fact that the 
d-wave order parameter changes sign at some regions in 
momentum space. Since in the d-wave case the order pa- 
rameter depends on momentum (or in other words angle 
9), a quasiparticlc will in general sense a different order 
parameter before and after scattering (normal scattering 
changes the momentum). When there is a sign differ- 
ence of the order parameter before and after scattering a 
zero-energy bound state appears (MGS). El These local- 
ized zero-energy states open up additional channels for 
current flow, leading to peaks in the conductance. Cal- 
culations of the current-voltage (I — V) cuxves for the 
normal metal-superconductor (NS) junctions confirmed 
the picture brought forward in Ref. [TI| MGS are in- 
finitely sharp when the transmission of the strength of 
the barrier is infinetely high. When the strength of the 
barrier decreases the MGS broaden. 

In addition, there are some other features related to 
the ac Josephson effect involving d-wave superconduc- 
tors that have not yet (to our knowledge) been ex- 
perimentally—established, such as negative differential 
conductancefijtJ and cancellation of the odd ac com- 
ponents, doubling the Josephson ffequencyJHI although, 
there are experiments discussing related features J3'E£l 
Negative differential concLuctance is a result of resonant 
conduction through MGSE3 Moreover, we will in this pa- 
per show an example of ZBCP for the case when there 
are MGS on both sides of the barrier (possible only when 
there are d-wave superconductors on both sides of the 
barrier). We note that for the cases we discuss ZBCP ap- 
pear for relatively small values of barrier strength. The 



geometry used in these calculations correspond to exper- 
iments just recently performed Jj'EEl 

Cancellation of the odd ac components is explained as 
follows (see Fig. ||). The keypoint is that the total cur- 
rent is a sum of contributions, associated with different 
angles 9. In the case of two d-wave superconductors with 
orientation angles a± = and oli = tt/A respectively, 
quasiparticles incident at angle —9 experience a different 
sign of the gap A 2 compared to quasiparticles incident 
at 9. This sign difference (which corresponds to a phase 
difference of it) enters the expression of ac components as 
e"™ 7r = (— l) m , where m is the index of the ac component 
in Eq. ([!]). This means that for even (odd) positive 
and negative angles add up (cancel out). 113 We note that 
the doubling of the Josephson frequency is completely in 
line with the change of the current-phase relationship in 
stationary Josephson junctions from being 23r-periodic to 
being 7r-periodic in some arrangements. l3~E3 

Moreover, there appears in voltage-biased junctions of 
d-wave superconductors a possibility to have current flow 
parallel to the junction. This is not necessarily related to 
the ac Josephson effect but could appear in NS junctions 
as well.113 




FIG. 2. Symmetry argument for the cancellation of the odd 
ac components for the specific orientation doMir/4- Consider- 
ing injection of electron-like quasiparticles at angles ±9, the 
gap in the left superconductor will be the same, while there 
will be a sign difference in the right superconductor. In the 
above figure, the negative sign appears at — 9 and can be 
thought of as an extra phase it for this injection angle. The 
current contributions from ±9 are equal in magnitude but 
have different signs for the odd components, because of the 
extra phase it, and will therefore cancel. 

Recently, there has been a discussion concerning the 
appearance of superconducting states breaking time- 
reversal symmetry, especially at inhomogeneities (like 
surfaces/interfaces) in d-wave superconductors. These 
states have been theoretically predicted to have an order 
parameter of the form d+is, where the zs-part is induced 
by the inhomogeneity and only present close to the in- 
homogeneity. The zs-wave part of the order parameter 
has vital physical importance: the d + is combination ex- 
plicitly breaks time-reversal symmetry, and peculiar sur- 
face/interface currents appear with the existence of the 
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is part. One example of this phenomenon is twin bound- 
aries of YBCO, whose order parameter is of the form 
d ± s (a real combination). The twin boundary divides 
the superconductor into two parts with different signs in 
front of the s-part. To match the two sides the most fa- 
vorable way (from a Ginzburg-Landau theoretical point 
of view) , close to the boundary a complex order parame- 
ter appears in order to avoid a discontinous change in the. 
phase between the s- and d-part of the order parameter 
This resutLhas also been found within a Green's function 
formalisnaea and by solving the Bogoliubov-de Gennes 
equations Another example is the surface of a rf-wave 
superconductor, where the surface lowers the symmetry 
in such a way that a combination of the subdominant s 
and the dominant d order parameters are now allowed 
at temperatures below a certain temperature defined by 
various properties, giving rise to a d+is surface state with 
associated time-revecsaJ-aymmetry breaking currents par- 
allel to the surface. Bc3~l3 These surface states are doubly 
degenerate, corresponding to two opposite current direc- 
tions. A third example (and most closely related to the 
case studied in this paper) is the junction between an s- 
wave superconductor and a d 2 .2_ y 2-wave superconductor 
rotated 45° ,E3 where a combination of s- and d- waves co- 
exists close to the junction area. Again the phase differ- 
ence of the energetically stable phase is such that a d±is 
order parameter appears, also this time doubly degener- 
ate corresponding to phase differences of ±7r/2 between 
the s- and ci-wave side. The same effect has theoretically 
been found in junctions between two d x 2_ y 2 supercon- 
ductors where one of the superconductors is rotated by 
45°. El 

All the examples discussed above concern the station- 
ary state. This is definitely different from our case w ith a, 
voltage present between the two sides of the junctionJiMZl 
In fact, when the junction is voltage-biased it does not 
care about the ground state phase difference between the 
two sides. Rather, all possible phases are taking part in 
an averaging procedure. Therefore, it is not surprising 
that the case discussed in Ref. |2t] (the third example in 
the previous paragraph) does not produce parallel cur- 
rent density for the zero-frequency ac component in the 
voltage-biased case, since the two degenerate states (cor- 
responding to two opposite current directions) are now 
averaged to produce a zero net dc current density. In 
Fig. |3| we show the mechanism in the voltage biased case 
for the net current density parallel to the interface be- 
tween the two sides. The bottomline is that a d x i_ y i- 
wave order parameter with orientation angle a^O treats 
quasiparticles associated with positive injection angles 
differently than quasiparticles associated with negative 
injection angles, leading to an imbalance between cur- 
rent contributions associated with positive and negative 
injection angles. 
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FIG. 3. This picture shows how the parallel (y- direct ion) 
current appears as a consequence of the asymmetry of the 
d-wave order parameter. On the left side we have a normal, 
an s-wave or a d-wave structure; on the right side we have a 
d-wave superconductor with a 7^ niv/A, where n is an inte- 
ger. Quasiparticles incident on the right side (from the left) 
in general experience a different gap depending on angle of in- 
cidence 9. Since current contributions (among other things) 
depend on the gap this means that the parallel current contri- 
bution corresponding to 9 does not necessarily cancel the one 
corresponding to —9. The solid and dashed lines illustrate 
electrons and holes propagating in both directions. 

Therefore, in this paper we consider the ac Joseph- 
son effect and focus on three main issues: first, current- 
voltage curves of d-wave S1IS2 junctions; second, the dis- 
appearance of the odd ac components; third, the current 
density parallel to the junction direction. Moreove r, . we , 
give a fairly complete account of the methods usedllMZl 
Especially, we elaborate on how various ways of writing 
down current formulas relate to each other. 

In the case of s-wave superconductivity Eq. ([!]) has 
been studied for a long time. In the tunneling limit 
(low transipission of the insulating barrier) detailed 
calculationstilLJ showed excellent agreement with exper- 
imental results for the conventional s-wave superconduc- 
tors. These calculations used the tunneling hamiltonian 
method. Due to the huge success of the tunneling calcu- 
lations less attention was paid to the non-tunneling case. 

The work by Blonder, Tinkham, and Klapwijk proved 
that with fairly simple methods it was possible to pro- 
duce current-voltage curves of NS junctions for any 
barrier strength.Lj Later, another line of research fo- 
cused on introducing—the concept of mesoscopics .iuJxt 
superconductivity.E^rEj Both thfvdr: Josephson effectEj'EEl 
and the ac Josephson effectLTEa were studied, again 
without the tunneling-limit restriction. 

The main tool in many of these papers dealing 
with the case of arbitrary transmission— has., beep the 
Bogoliubov-de Gennes (BdG) equation.EJQE3'El The 
time-dependent BdG equation together with a formula 
for the current gives a complete description of the cur- 
rent in both the stationary and non-stationary case. In 
this paper we use a generalizatipnliaO of the methods 
worked out for the s-wave case.E3cil 

When an electron with subgap energy from a normal 
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region approaches the superconductor, there seems to be 
no way of transmitting electric charge into the supercon- 
ductor: there are no propagating states of the supercon- 
ductor in the subgap region. This becomes evident when 
solving the Bogoliubov-dc Gcnncs (BdG) equation result- 
ing in d ecayi ng wave functions when solving for subgap 
energies E5e3 

However, there is one important physical process mak- 
ing transmission of electric charge possihle_also for sub- 
gap incidence: Andreev reflection (AR).oEj In this pro- 
cess an electron (from a normal region) with subgap en- 
ergy incident on a superconductor is reflected as a hole 
moving in the opposite direction to the original electron. 
This event does obviously not in itself conserve particles. 
Therefore, a Cooper pair of electrons is injected into the 
superconductor. When two superconductors embrace the 
normal region, forming an SNS junction, the hole will fi- 
nally hit the other superconductor and again Andreev 
reflects, this time as an electron. These AR's will con- 
tinue back and forth in the normal region. Since the 
voltage drop is over the normal region, the electrons and 
holes gain the energy eV every time they pass the normal 
region (see Fig. ||). These repeated reflections are usually 
named multiple Andreev reflections (MAR). The reflec- 
tions are perfect inside the gap but decrease in strength 
when the magnitude of the excitation energy is increased 
away from the gap region. 

For each round trip two AR's transform the original 
electron into an electron at an energy 2eV higher than 
from the start. This means that the wave functions con- 
tain parts (sidebands) separated by multiples of 2eV from 
each other, leading to an ac current with a frequency of 
2eV/h as seen in Eq. (|l|). 

We think Fig. | is a key figure in order to illustrate 
the underlying physical mechanism of the phenomenon 
dicussed in this paper (the ac Josephson effect): it makes 
clear the meaning of Eq. (Q). There are two scattering 
events of our problem: first, normal scattering at the 
barrier described by a transmission amplitude; second, 
AR at the NS interfaces described by an AR amplitude. 
Therefore, Eq. (|l|) can be thought of as an expansion 
in transmission amplitude (through the barrier) and AR 
amplitude. Higher orders of the wave function corre- 
sponding to sidebands visualized in Fig. |^ contain higher 
orders of transmission through the barrier and higher or- 
ders of AR. 
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FIG. 4. The scattering state originating from an elec- 
tron-like quasiparticle incident at energy E. There are elec- 
trons and holes propagating in both directions (no arrows 
included in the figure) in the normal region at energies 
E„ = E + neV, where n is an integer and V is the volt- 
age over the junction. The shaded area is the junction region 
where the voltage drop occurs. 

The present paper is organized in the following way. 
In Section [n] we outline the model, which is used in Sec- 
tion III to calculate the various ac components of the 
We present the results in Section IV, and con- 



current. 

elude with a summary in Section |v} We feel that an im- 
portant part of our paper is given in the appendices: in 
Appendix |A] the time-dependent BdG equation is solved, 
and we also give some numerical advice on how to calcu- 
late the current; in Appendix [B] we compare how various 
current formulas previously used in the literature relate 
to each other; in Appendix ^| we point out a couple of 
symmetry relations that are relevant to our problem; fi- 
nally, in Appendix [5] we write down a pair of analytical 
expressions describing the ballistic case. 



II. THE MODEL 

To mathematically describe the voltage-biased Si | S2 
system shown in Fig. [I] we have made a number of as- 
sumptions, some of them pertaining to the superconduct- 
ing parts of the system and some of them pertaining to 
the normal part of the system. In the following we detail 
these assumptions. 

Starting with the normal part (the junction) in be- 
tween the superconductors, we assume a planar struc- 
ture with translational invariance in the y-direction. The 
planar geometry is a necessary condition for_MGS to ex- 
ist, contrary to the case of point contactsa where no 
MGS are around. Moreover, we neglect in our treat- 
ment any influence from interface roughness. In the nor- 
mal region we only take into account scattering due to 



4 



the barrier, which is modeled as a one-particle potential 
V(x) = H5{x). The 5-function form of V(x) results in a 
reflection amplitude r(6) = Z/(icos8 — Z) and a trans- 
mission amplitude t{9) = Laos 6 /(i cos — Z), where Z is 
defined as Z = 2mH/h 2 .c3 The parameter Z describes 
the strength of the barrier: Z = corresponds to no 
scattering at all (ballistic case); Z = oo corresponds to 
the case when the two sides of the junction are uncoupled 
(no x-current in this case); a high value of Z corresponds 
to the tunneling case. 

Next we consider the superconducting parts. In our 
calculation self-consistency is not fulfilled. In the liter- 
ature there is work that has included self-consistency in 
the tunneling limit showing, that non-zero energy states 
appear along with MGS.Ii-3'EjIn our work we neglect these 
effects; on the other hand, we allow the barrier strength 
Z of the barrier to have any value. 

Therefore, we assume that the superconducting order 
parameter of the structure in Fig. |l| is 

A _fA 1 (fl), x<0 

A ~ \ A 2 (0)e^°, a;>0. (Z) 

Since the overall phase is unimportant we can here choose 
Ai real. There may be a phase-difference </>o over the 
junction which we defer to the right-hand superconduc- 
tor. The angle 8 corresponds to the angle of incidence of 
the quasiparticle approaching the superconductor. In the 
following we will use various forms of order parameters 
for Ai(9) and A 2 (#) ranging from the s-wave case to the 
d-wave case. The d-wave superconductors are allowed to 
be rotated an angle a measured relative to the interface 
normal as described in Fig. |l|. 

In this paper we do not include the magnetic field. In- 
cluding the magnetic field will screen any current on the 
length scale of the London penetration depth, forcing the 
current to flow next to surfaces/interfaces of supercon- 
ductors. Neglecting screening effects is traditionally con- 
sidered valid if the system under study has dimensions 
smaller than the London penetration depth. However, a 
finite system immediately raises questions about where 
the parallel current finally ends up. In this paper we do 
not address this question; instead we just note that a full 
treatment including screening effects will be required to 
understand how the parallel current is guided close to the 
interface of the d-wave superconductor. 

Using the models discussed above for both the nor- 
mal and the superconducting regions we solve the time- 
dependent Bogoliubov-de Gennes (BdG) equation for 
anisotropic superconductors 




/ ho(x,x') A(x,x',i) \ /u(x',t)\ _ 
^ A*(x,x',i) -/i (x,x') J {v(x',t) J 



m (t$)> Vx,x')=*(x-xO(f|-;), (3) 

where fi is the chemical potential, and A(x, x',i) is the 
order parameter. Introducing a voltage the chemical po- 
tentials of the two reservoirs will be shifted relative to 



each other, eV = fi2 — The natural energy refer- 
ence level in a superconductor is the chemical potential. 
In order to establish a global energy reference level, a 
gauge transformation is performed as described in detail 
in Ref. [li| As a result of this transformation a phase dif- 
ference appears over the junction obeying the Josephson 
relation dt<fr — 2eV/%. In other words, for our case with 
constant voltage we describe this time-dependent phase 
difference multiplying A 2 (#) by exp(i2eVt/?i). 

In Appendix ^we solve Eq. (|^) for the wave function 
= (V'vjVu) corresponding to an electron- like quasi- 
particle incident on the junction at an energy E v . Using 
this wave function we will in the next section calculate 
the current. 



III. CALCULATING THE CURRENT 

To calculate the current density j we use the general 
formula 

j = — V [2f(E u ) - 1] Im {u* v Vu u + v* v Vv v } , (4) 

771 *■ — ' 

v 

where we introduced the Fermi distribution function 
f(E) = 1/[1 + exp(E/k B T)}. The sum is over all scat- 
tering states originating from electron-like quasiparticles 
approaching the junction from the superconducting reser- 
voirs at both negative and positive energies E v . There 
are different ways of writing the formula for the current 
density, but as shown in Appendix [b] they are all equiv- 
alent. Note that j in our case may have a component 
in the y-direction (parallel to the junction). The wave 
function coefficients determined in Appendix ^ are in- 
troduced into the current formula of Eq. (^) in a way 
outlined below. 



A. The a;-current 

We start by studying the current density j x perpen- 
dicular to the junction. Since we consider a structure 
smaller than the London penetration depth (allowing us 
to neglect self-field effects) the current density will dis- 
tribute uniformly over the junction. Therefore we may 
calculate the current density at any y-coordinate. More- 
over, a conservation law for the current density assures 
that we may calculate j x at any ^-coordinate; it proves 
convenient to calculate the current density in the nor- 
mal region (in our treatment to the left of the barrier 
in Fig. [l]) between the superconducting reservoirs. The 
total perpendicular current I x per a6-plane is then given 
by Ix = Lyjx, where L y is the length of the junction in 
the j/-direction. 

In Eq. (^J) the quantum number v is the wave vector k 
of an electron-like quasiparticle approaching the junction. 
The direction could be thought of as another quantum 
number; we introduce r = + (— ) to indicate that the 
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quasiparticle is incident from the left (right). Therefore, 
using the wave function coefficients of Appendix ^[ the 
current per 06-plane is 

4 = Ly — -—^-— ^ TCOs(# k )T T (k)[2/(£' k ) - 1], (5) 
m Li x -u v ' — 

where the combination olkp (the Fermi wavevector) and 
cos(6*k) (6*k is the angle of incidence for an electron-like 
quasiparticle incident with momentum k) is produced by 
the x-derivative in Eq. (Q). The factor l/(L x L y ) in the 
equation above normalizes the wavefunction of the in- 
coming electronlike quasiparticle. In Eq. (^) we have 
introduced 

T T (k) = [(aE,2n+2m)* a L,2n ~ ( rf L,2n+2m)* rf E,2ri 

~K&£,2n+2j7i) ^L,2n — ( c L,2n+2m) C L,2n] X 



^irm2eVt / h ^irm^o 



(6) 



The set of states (electron-like quasiparticles) to sum over 
in Eq. (Q) is illustrated in Fig. [ll](c). In Eq. (^]) n and m 
are integer numbers. 

Next we change the sum over k in Eq. Q into an inte- 
gration over d 2 k. This produces a factor of L x L y /(2Tr) 2 . 
Finally, changing the integration over d 2 k into an integra- 
tion over the superconducting energy E = [(h 2 k 2 /2m — 
/i) 2 + A(6' k ) 2 ] 1 / 2 and the angle of incidence 9 = 6> k , we 
find the following expression for the current 



ix_ev j_ r /2 

V Q ~ + 4D J_ w/ 



/it/2 i-OO ri- 

d9cosfl / — [2/(i?)-l] 
-tt/2 J-oo A 



O"0 — Ly 



2S/ 2 em 1 / 2 E 1 F /2 A D 



h 2 



(7) 



where D = j d9\t(9)\ 2 cos 9/2 and is the Fermi en- 
ergy. The first term proportional to eV in Eq. (0) is 
due to refering energies of right-movers (left-movers) to 
the chemical potential of the left (right) superconductors. 
For details, see Refs. [l4] and ^6|. We note that the sec- 
ond term of Eq. (Jfy is zero when all order parameters are 
put to zero. It is therefore reasonable to call this part 
the superconducting contribution to the current and the 
first part the normal contribution to the current. 

It is now convenient to introduce f2 m = m2eVt/h + 
m(j)Q. We then rewrite Eq. (Q) 



cr A ^ 



cos fL 



B x ,m sin r2 m ] , 



(8) 



where A x>m and B X Tn are given by 

i r' 2 

A x , m (y) = — / d9 co 

iLJ J -tt/2 



f°° dE . 


~ E \ 


/ -— tanh 




l-oo A 1 


K 2k B T) 



x^ T N T (E)Re{T^E,9,m)}, 



1 /*7T / 2 /» OO / ^ 

B„„(>')--^| i/2 <i<'» 8S y_-ta„ ll ( — 

x^« r (E)Im{r;(E,8,m)}, (9) 

r 

together with 

OO 

T£(E,9,m)= ^ [( tt 2n+2m)* ffl L " (4n+2m)* d 2ii 

TL— — OO 

+ (^2n+2m) ^2n — ( C 2n+2m) C 2n 



*cLI. (io) 



In Eq. (||) we write 2/(^7) - 1 = tanh(-£/2fc B T). 



Final 
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y, we introduce an amplitude C x<m = (A 2 

Bx.m) 1 ^ 2 an d a phase a 1:m = arctan(i? 2 , : „ l /A 2 , : „ l ) in or- 
der to write the current as 



* \ / ^x,m COSy\L m OL x ,m) 
An 



CO 



O i(flm-o, !m ) 



(11) 



where we in the last step used the symmetries A x ^ m = 

Ax,mi ^x,—m ^4,m giving C x — m — C x ,mi &x,—m 

—OL Xi ra. These symmetries are due to the relation 
T T (E,9,-m) = \T T (E,e,m)}*. Note that the last ex- 
pression in Eq. ( |lT| ) is equal to Eq. ([[]). 

Studying the phases Vt m — a Xtm in Eq. ([ll]) we see 
that for to = they are all zero implying that the 
zero-frequency component does not depend on the phase 
difference (f>o over the junction. For the components 
with m ^ 0, we see that we can arrange the phases as 
^ m -a x . m = mojj(t-t ), where t = (a x , m /m-(f> )/ujj 
can be thought of as the time when the measurement 
starts. It will therefore be of no importance for the differ- 
ent current components in the following discussion, which 
rather focuses on the amplitude C m . 



B. The ^-current 

We now move on to study the y-current (current paral- 
lel to the junction). Below we will see that in general the 
y-current density j y (x) (current per afo-plane and unit 
length in the x-direction) depends on the x-direction. To 
stress this dependence on x we have chosen to discuss the 
parallel current in terms of the current density j y . 

The analysis of the parallel current density follows 
closely the one presented for the perpendicular current 
above. There are four main differences: first, the deriva- 
tive with respect to y in Eq. (|J) now leads to a factor 
of sin# (instead of cos 9); second, there will only be plus 
signs in the equation corresponding to Eq. (|l0|); third, 
there is an oscillating term in j y ; fourth, there is no nor- 
mal contribution to the parallel current as in the perpen- 
dicular case discussed above. 
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Therefore, we write down (suppressing the IE- 
dependence) the parallel current density j y as 



Jy 



imujjt 



(12) 



with the components defined as 

JvAv) 



AyAV), 

C v , m {V)e l ^'- a y^\ m^O 



((T0/Ly) 

C v , m (V) = yJ^, m (V) + Bl m (V), 

By y m(V) 



a y,m — arctan , 



(13) 



together with 

x^TW r (E)Re{IJ(E,9,m)}, 

T 

B y , m (y) = — — / d6> sin 6> / — tanh 



AD 



Tl/2 



A 



2k B T 



xJ2N T (E)lm{T^(E 7 9,m)} , 

T 

OO 

r;(£, 9,m) = J2 {(«2„ +2m )*«2„ + (^2m)% 



71— — OO 



+ ( c 2n+2m) c 2n + (^2ri+2m) ^2rt 

+ 2 [( a L,n+m)*^L,ri + (^.ri+m)* !,™] 
2i/cj^ cos f?x ^ 



xe 



(14) 



The last term in Ty(E,6,m) oscillates on a scale 1/kp 
and would be hard to observe in an experiment. We 
therefore drop it when we calculate the parallel current. 
This is reasonable since j y is a current density: To get the 
total current we should integrate in the x-direction. The 
oscillating term will then average to zero. Dropping the 
oscillating term, j y will no longer be continuous when 
passing the barrier: the value of j y will be different at 
x = + and x = 0~ (with the (^-function barrier located 
at x = 0). However, we have checked that the m = 
component of j y is continuous if the oscillating term is 
included. For completeness we will in this paper give 
(the average value of) the current density on both sides 
of the barrier, denoted j y L and j y n, respectively. For the 
ballistic case (Z = 0) we have jyL — JyR- 

In the next section we discuss C x<m (V) and C y<m (V) 
for some configurations of superconductors. 



IV. RESULTS AND DISCUSSION 

We will in this section present curves only for positive 
voltage: in Appendix |c| we show that the current is odd 
in voltage. 

We start by discussing the results for the m = term 
of the sum in Eq. (|l|) . In Fig. || we have plotted the I — V 
curves for different values of the temperature. For zero 
temperature there is a peak in the current for some-Otis 
entations due to resonant transport through MGS.EjEj 
This feature (resulting in negative differential conduc- 
tance) is especially pronounced for arrangements where 
one superconductor is an s-wave superconductor and the 
other one is a d-wave superconductor with orientation 
angle a — tt/4. In Fig. |](a)-(b) we study how temper- 
ature influences the current peak appearing at a voltage 
eV = A s . We find that the peak is rather robust up 
to relatively high temperatures, such as 0.5T CS shown in 
the figure, where T cs is the critical temperature of the 
s-wave superconductor. In Fig. || we have taken into ac- 
count the fact that the s-wave superconducting gap is 
reduced when the temperature is increased according to 
the theory by Bardeen, Cooper, and Schrieffer.cZl This is 
the reason for the shift of the peak towards lower voltage 
when the temperature is increased. The d-wave order 
parameter, on the other hand, is assumed not to be af- 
fected by the temperature variations, since T c for HTS 
compounds in general is considerably higher compared 
to conventional superconductors. ■— . 

It is interesting to compare with the NS junction,E3 
where it has been found that the presence of MGS induces 
a dramatic increase of current for low voltage resulting 
in ZBCP. This structure is more sensitive to temperature 
than the current peak of the s|d^/4-junction discussed in 
the paragraph above. The reason for this greater sensi- 
tivity is the density of states: it is well known that MGS 
are located at the Fermi energy on the (i-wave side. This 
means that quasiparticles incident from the normal side 
will reach MGS for any small voltage. This, however, is 
in great contrast to the s|d T /4-case, where MGS on the 
d-wave side may only be reached from the s-wave side 
at a voltage greater than A s simply because the super- 
conducting density of states is zero for subgap energies. 
The conclusion is therefore that the NS case will be more 
affected by thermal smearing introduced by the Fermi 
functions, since this smearing is centralized around the 
Fermi energy. This argument breaks down at tempera- 
tures larger than the s-wave gap, taking place close to T c 
of the s-wave superconductor. 

An interesting paperQ has recently appeared studying 
experimentally the system corresponding to the case de- 
picted in Fig. |^. In that paper two sets of curves are 
presented: first, Fig. 4(a) of Ref. denoted here by (i); 
second, Fig. 4(b) of Ref. R denoted here by (ii). While the 
results of (i) are in agreement with the picture outlined 
in our Fig. H, the curves of (ii) are not easily understood 
from our calculations. The critical point for not having 
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agreement is the fact that in (ii) there is no opening of 
the superconducting s-wave gap corresponding to the Pb 
electrode to be found. 
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FIG. 5. Current-voltage curves and differential conduc- 
tance-voltage curves for s\d^/4 [(a) and (b) respectively] and 
N\d 7I /4 junctions [(c) and (d) respectively] are shown for some 
temperatures. The zero bias conductance peak present for 
the N\d case is moved to finite bias eV = A S (T) when the 
counter electrode is an s-wave superconductor. In addition 
negative conductance appears. The peak moves to smaller 
voltage when T is close to T cs since A s is reduced by temper- 
ature according to the BCS theory. The figure demonstrates 
that the peak is less sensitive to temperature for the s\d junc- 
tion compared to the N\d junction. The barrier strength is 
Z=5, and A S (T = 0) = 0.2A . 



Very recently an investigation of ZBCJw-jp bicrystal 
grain boundary junctions have appeared,BH3 discussing 
(among other things) a geometry corresponding to a.\ = 
— «2 = 18.4°, in which ZBCP was found. In our calcu- 
lations ZBCP can only be found in rather transparent 
(small Z) d\d junctions, as shown in Fig [| The reason 
for this is that the scattering states start from the bulk 
where there are no superconducting density of states close 
to zero energy except for angles corresponding to the gap 
nodes of the ci-wave superconductor. For small voltage, 
the MGS can only be reached by the scattering states 
from those angles, making the MGS contribution very 
limited in the angle average. Broadening of MGS due to 
reduced barrier height Z (which is incorporated in our 
treatment from the start) increases the range of angles 
from which the scattering states may reach MGS. For low 
barrier height the width of the MGS becomes, comparable 
to the gap making the junction ballistictZlc3 for scatter- 
ing states starting near the gap edges, allowing MAR to 
contribute massively to the current. The current then 
takes a finite value for small voltage and a huge ZBCP 
appears. It is important to have MGS on both sides of the 
barrier for this effect; otherwise, the number of AR will 
be reduced, resulting in a smaller current contribution. 
Fig. H demonstrates that at a barrier strength Z = 1, cor- 
responding to an averaged transparancy D = 0.38, the 



ZBCP has disappeared within our model: for this barrier 
strength the width of the MGS is not high enough to give 
large current at small voltage. 
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6. Current- voltage curves are plotted for the 
n /s junction with barrier strenths Z — 0, 0.25, 0.5, 
and 1. The finite value of the current at small voltage for 
low (but finite) barrier strength is due to resonant transport 
through MGS making the junction ballistic. For comparison, 
the ballistic (Z — 0) curve is shown (no MGS in the ballisic 
case) where there is always a finite current in the zero-bias 
limit. 

A related issue is the case ot\ = a 2 = 45° discussed 
in Ref. 13. Also in this case a combination of broadened 



MGS on both sides of the barrier and MAR give rise to 
ZBCP for small (e. g. Z = 0.25) barrier strengths, al- 
though the peak is not as sharp as for d n /g,\ d- n /g,. As 
explained above there are no scattering states at zero 
energy except for the angles corresponding to the gap 
nodes. Since for the orientation ct\ = a 2 = 45° the gap 
nodes are geometrically lined up against each other, the 
influence of MGS will be less dramatic compared to the 
case discussed in the previous paragraph. An interest- 
ing paper, mainly on calculations in the tunneling limit, 
introduces phenomenologically an imaginary part pf the 
energy, corresponding to broadening of the MGS.lj (In 
our calculation there is no imaginary part of the energy.) 
Due to the imaginary part ZBCP was found in the case 
Q'i = a?2 = 45°. E3 The difference compared to our work is 
that in our work the broadening is not introduced phe- 
nomenologically; contrary, it is microscopically described 
in terms of reduced strength of the barrier. 

We now continue to discuss the m ^ ac current com- 
ponents of Eq. (Q). In Fig. |?] we show some results for 
current density for both the perpendicular (x) and par- 
allel (y) direction. The curves shown in Fig. |?] are for 
the ballistic [Z = 0) case. The ballistic case is rather 
restricted; however, it will prove very useful for us to 
illustrate some important properties of d ai= o\d a2 junc- 
tions in the ac case. These qualitative properties carry 
over to the non-ballistic case in a straightforward man- 
ner. Quantitatively the curves will look different in the 
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presence of normal scattering. 

We start by discussing ct2 — 7r/4. In Fig. we plot the 
magnitude C m of the ac components. It is then possible 
for us to show how the -disappearance of the ac compo- 
nents in the x-directionli3 will be complemented by the 
appearance of the corresponding ac components in the 
y-direction as seen in Fig. 0(b) and (e). The detailed 
explanation for this effect is given in Appendix |c| and 
is good for any value of Z. The main point is that the 
current is calculated from an integration over angles 9, 
where the integrand is composed of two functions. The 
first function is cos(0) [sin(#)] produced by a derivative 
in x (y) in Eq. (Q). The other function can be proved to 
be odd (even) in 9 for odd (even) m. It is now clear that 
the combination of these two functions integrated will re- 
sult in zero current density in the x (y) direction for odd 
(even) m. This scenario is confirmed in detail from the 
calculation presented in Fig. ^. We note in passing that 
the m = component has a finite current in the limit 



eV 



as discussed in Appendix jp| 
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FIG. 7. The first three components (m — 0, 1, 2) of the per- 
pendicular current [(a)-(c)] and the parallel current [(d)-(e)] 
is shown for the do\d a2 junction for three different orienta- 
tions of the right superconductor. The odd components of 
the perpendicular current is suppressed when the orientation 
angle oli is changed from towards 7r/4 (and zero at 7r/4). 
At the same time the odd component of the parallel current 
is enhanced. We also see that all components of the paral- 
lel current vanishes for the orientation ct2 = 0. The barrier 
strength is Z = and the temperature is zero. 

Next we discuss the parallel current density in the pres- 
ence of normal scattering. We specialize to the do |cL/8 
junction with Z = 1 for zero temperature. The curves 
depicted in Fig. ^ are the m = component calc ulated in 
the normal regions. As explained in Section III the aver- 
age current density in the left side differs from the right 
side of the barrier. We note that the large-voltage limit 
of the curves in Fig. is of the order of the perpendicu- 
lar excess current density, implying a saturation at large 



voltages. Below we indicate which side the current den- 
sities are to be associated with by adding an index L/R 
for left/right whenever the side of the junction matters. 

The fine structure in Fig. || is due to subharmonic gap 
structure of the-same origin as in the case of perpendic- 
ular current .E3EZI 
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FIG. 8. Parallel current for the do\d^/ s junction to the left 
and to the right of the barrier. The barrier strength is Z=l 
and the temperature is zero. 

The large-voltage limit of the parallel current density 
can be calculated analytically in the ballistic limit, see 
Appendix [d| For the orientation ai = and a.2 = 7r/8, 
the value of the large- voltage limit of the parallel current 
density will be negative. In Fig. |we plot the dependence 
of the large-voltage limit of the parallel current density 
on the barrier strength Z . The result is that we get a sign 
change at a rather small value of Z. The sign (change) 
can be qualitatively understood by studying Fig. |3| in 
combination with Eq. (Dl). We see that Fig. ^ approx- 
imately represent the case = tt/8. For Z = the 
current contributions for each angle of incidence will be 
averaged by the function sm9. We see in Fig. || that 
the largest gap appears at negative angles where sin 9 
is negative. The value of the integral in Eq. (|Dl]) will 
therefore be negativ e, sin ce only the absolute value of the 
gap appears in Eq. (Dl). Introducing normal scattering 
means that small angles of incidence are more strongly 
weighted (note the structure of transmission amplitude 
t(9)). In Fig. H we see that the largest gap will be at pos- 
itive angles in this case, implying a positive current since 
sin 9 is positive for positive angles. In addition, the nor- 
mal scattering will introduce contributions from the gap 
A2 = A2(— a-i) at the expense of A2. One can show that 
for the do\d a2 junction the relation j y (— a-i) — —jy(&2) 
holds for all Z, implying that the introduction of A2 
gives contributions with the positive sign. The inset in 
Fig. ^illustrates the relation j y {— 0.2) = — jyip-2)- In con- 
clusion, the positive contributions will quickly dominate 
over negative contributions in the angle average when we 
change Z from zero to finite values. 
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FIG. 9. We show the current density parallel to the junc- 
tion interface in the limit eV >> Ao as a function of the 
barrier strength for the orientation «i = and Q2 = 7r/8. 
In the inset we show the dependence on the misorientation 
angle «2 of the right superconductor for the barrier strengths 
Z — and Z = 1 keeping ot\ = 0. Note that the factor D 
(which is Z-dependent) is included in ao, giving a larger nor- 
malized value of the current for non-zero Z. Zero temperature 
is assumed in all cases. 

Finally, we show in Fig. [l0| some I — V curves for the 
N\d a junction for several barrier strengths and orienta- 
tions a. The curves demonstrate what is illustrated in 
Fig. U for almost any orientation (the exceptions are 
a2 = twt/4) there is a net parallel current density. In our 
model the parallel current density will be present every- 
where in the normal region. This might not be surprising 
since we do not take into account any impurity scattering 
in the bulk. 
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FIG. 10. Currents for the N\d a junction at barrier 
strengths Z — and Z = 1. In (a)-(b) we show the perpen- 
dicular current I x flowing through the junction and in (c)-(d) 
we show the parallel current density j y . In (b) we see that 
the zero bias conductance peak grows up when we rotate the 
orientation of the superconductor from Q2 = to Q2 = 7r/4. 
The parallel current only exists when 02 7^ mr/4, as seen in 
(c)-(d). 



An important limitation of the present work is that we 
do not include magnetic field. From a calculation with a 
self-consistent magnetic field one should expect that the 
parallel current density is screened out in the supercon- 
ductors on the length scale of the London penetration 
depth. This means that our results for the parallel cur- 
rent density presented in Figs. |?] to |l0| are only valid for 
systems with dimensions smaller than the London pene- 
tration depth. 



V. SUMMARY 

In this paper we have studied the ac Josephson effect 
for planar superconducting d-wave junctions. The pres- 
ence of zero-energy midgap states (MGS) at interfaces of 
e?-wave superconductors for some orientations influence 
the current-voltage curves (the zeroth order ac compo- 
nent), since MGS open up additional channels for cur- 
rent flow. We have primarily discussed two aspects of 
MGS in this paper. First, our results feature current 
peaks at finite voltage (resulting, in negative differential 
conductance) in the Si | S2 case.113 We have in this paper 
explored how sensitive this current peak is to tempera- 
ture. We have also discussed how well our calculations 
compare with an experiment :I3 some curves of this exper- 
iment compare well with our results; however, there is in 
Ref . a set of curves which we find difficult to understand 
within the framework of the model we are using. The 
second aspect relates to a recent experiment on d a \d- a 
j unctions £rE3 where zero-bias conductance peaks (corre- 
sponding to large current increase for small valtage) have 
been found. Our model explains these experimental re- 
sults in terms of resonant transport through MGS if the 
grain boundary is rather transparent, allowing multiple 
Andreev reflections to contribute hugely to the current 
for low voltage. 

An interesting feature of superconducting d-wave junc- 
tions is the possiblity to have parallel current density at 
the interface of the two electrodes. We explain the reason 
for this effect for the ac Josephson effect. We also demon- 
strate how the disappearance p£ the odd ac components 
for the perpendicular currentEZI in some orientations is 
complemented by appearance of odd ac components in 
the direction parallel to the junction. 

Moreover, we have given a fairly complete account of 
our method which is applicable to all types of junctions: 
all the way from the ballistic (no barrier) to the insulat- 
ing case. We have also discussed how different current 
formulas previously used relate to each other. 
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APPENDIX A: SOLUTION OF THE 
BOGOLIUBOV-DE GENNES EQUATION 

The complicated integral operator makes exact solu- 
tions of Eq. (||) not easily accessible. Therefore, we solve 
Eq. (||) approximately by an ansatz in each region of the 
following type: 



(v(x,t) J -Z. Ake U J 



-iE k t/h 



(Al) 



where the Fourier sum over different frequencies reflects 
the fact that there is a time dependence in the order 
parameter after a global reference level is chosen as de- 
scribed in Section [n} We follow the standard treatment 
when approximating the integral operator of Eq. (|j|): the 
center of mass coordinate R = (x + x')/2 and the rel- 
ative mass coordinate r = x — x' are introduced. This 
means that in these new coordinates the order parameter 
is written A(x, x') = A(r, R). The integral operator of 
Eq. (|J) is now approximated as 

/^A(x, x 'Mx-, ( ) = / d rA(r,x-r/2 W x-M) 
« £ A(k, x)A kWk e lk - x e- ljEk * /?i , (A2) 

k 

where we have introduced the Fourier transform A(k, R) 
of the order parameter A(r, R). In Eq. ( A2) we approx- 
imated A(r,x — r/2) w A(r,x), valid since A(r,x) is 
assumed to be a slow (fast) function in x (r). Since in 
the weak-coupling limit the pair potential is expected to 
be non-zero only close to kp , one can replace the momen- 
tum by an angle 9. The approximation made in Eq. ( A2) 
neglects terms of the order (Jcf£o) _1 and is called the 
quasi-classical approximations 

Before solving the time-dependent BdG equation it is 
convenient to introduce some notational simplifications. 
First we define the BCS coherence factors 



Vi(E,e) = E-sgn(E)^/E 2 -A l (9)' 
Ui (E,6) Ai(6) 
Vi(E,e) _E-i^/A l (9) 2 -E 2 



, \E\ > \M8)\ 



Ui(E,( 



A;(0) 



, \E\ < \Ai(9)\ (A3) 



where \ui\ 2 + \vi\ 2 — 1 and i — 1,2 refers to supercon- 
ductor 1 and 2. Using this definition, one can express 
the Andreev reflection amplitude at superconductor i as 
Ai(E n ,9) = A %n = v i (E n ,9)/u i {E n ,9), where E n = 
E+neV. The transmission amplitude for an electron-like 
quasiparticle at energy E and angle 9 from superconduc- 
tor i to enter the electron branch in the normal region is 



J t {E,B) = Ji = {u 2 (E,9)-vf(E,9))/ Ul (E,9). Animpor- 
tant notational detail in this paper are the barred quan- 
tities: these are to be associated with the angle 9 = ir — 9. 

To solve the time-dependent BdG equatian-wcJbllow 
the method worked out for the s-wave cascrffit^rrlcJ ; 



and 

further generalized to the anisotropic case£3 ! EJ We use 
the boundary condition that far away from the junction 
the wave function corresponds to an electron-like quasi- 
particle incident on the junction from one side at an en- 
ergy E and an angle 9. Thanks to symmetries of the 
BdG equation it is enough to study these electron-like 
quasiparticle solutions as explained in Appendix |b[ The 
incident quasiparticle gives in general rise to both re- 
flected and transmitted waves. Therefore, we start writ- 
ing down the wave function in the left superconductor as 
(corresponding to a quasiparticle incident on the junction 
from the left) 



^ ol v 1 (E n ,9) 1 ' 



I Ul (E n ,9) 1 



zq x 



di 



ui{E n ,6) 
vi(E n ,0) 



Bnt,n±0 n 



(A4) 



where q e /' 1 corresponds to the wave vector of 
electron/hole-like quasiparticles. For the wave vec- 
tors above (and in the following) we use the notation 
q = _(q x ,q.y) = q(cqs9,shi9) and q = {-q x ,q y ) = 
q(cos9,sm9), where 9 = it — 9 corresponds to the an- 
gle after normal scattering. We will approximate the 
magnitude of any wave vector in this paper as the Fermi 
wave vector, q « kp. Again, all barred quantities are 
to be associated with the angle 9. In the voltage case 
the injected quasiparticle goes through multiple Andreev 
reflections (MAR). The wavefunction therefore contains 
sidebands at energies E n = E + neV, where n is an even 
integer. 

In the normal region the wave function to the left and 
right of the barrier is 



E 



zk e -x 



aL, n e lk x + d L ^ n e _ 
bL, n e* h - x + c L:n e* h -* 



, (A5) 



^-t _ I a R ^ n e lW ' x + d R . n e 



ik e -x 



-*(-r+-§ tt ) ) ( A6 ) 



respectively. The wave function of the right supercon- 
ductor is 



So 



E 



a 2 , 



u 2 {E n ,9)e^e^) 



v 2 (E n ,9)e- 



_i£o _i(eVt) 



C2 ,■ 



v 2 (E n ,9)e l -e<- 



U2(E n ,9)e~ 



xe 



(A7) 
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In the wave functions written down in Eqs. (A4)-(A7) 
we have explicitly factored out the dependence on the 
phase 0o for convenience, making the coefficients a, b, c, 
and d all independent of 4>q. In the final expression for 
the current, however, the phase </>o will reappear. 

In the following we solve the matching problem. We 
have three interfaces where the wave functions and 



derivatives in Eqs. (A4)-(A7) are to be matched. The 
procedure is the same as the one carried out in Refs. ^ 
and [l6]. When matching at the left interface between 
superconductor 1 and the left normal region we get 



O-L.n — ^l,nbL,n + ^l#n,0) c L,n 



(A8) 



after eliminating coefficients pertaining to the left super- 
conductor. 

Repeating the procedure at the interface between the 
right normal region and superconductor 2 we have 



5R,ti+2 



A 



2,n+lCR,n+2 



(A9) 



We have now eliminated coefficients related to the super- 
conducting regions. 

It now remains to match the wave functions in the 
normal region at the scattering center. This procedure 
may be expressed in terms of the scattering matrix 



S e (0) 



r(8) t(8) 

t{6) -r*(6)t(6)/t*(0) 



(A10) 



The scattering matrix for holes is Sh(0) — S*(9). The 
scattering matrix relates incoming waves with outgoing 
waves. As outlined in Ref. [14] for the case of s-wave 
superconductors with different gaps, the wave-function 
coefficients of the left and right side of the normal region 
fulfill [suppressing the 0-dependence of r(9) and t(9)] 



dh,n = rAx.nbL^n + rJi<5„, + tA 2t n+lCRM+2, (All) 

r*t - 

0,R,n = tAi^ n bL,n + tJlS n fi — A 2 ,n+lCfl iJl+2 , (A12) 



t* 



rt* 



CR.n = t*Ai n dL,: 



t 



■ A-2 ,Ti-lflfl,n-2- 



(A13) 
(A14) 



Using Eqs. (A8)-(A14) we solve for d^. As a result we 
find the following recursive relation: 



a n dL,n+2 + Pnd-L ,71 i Tn^L,n— 2 

where 

^2, 71+1^1, 71+2 



rJlSnA 



(A15) 



a n = - t 



1 — ^2,71+1^4-2,71 + 1 

J / ^1, 71^1, 71 



f3 n = l- A hn A hn + \t\ 

^-2, 71+1^-2, 71+1 



1 — ^2,71-1^42,77,-1 

^4.2,77-1^4.1,77, 



(A16) 



7» = -t 



1 — ^-2, 71 + 1^-2, 71 + 1 / 1 — ^2, 71-1^-2, 71-1 

The coefficient is determined from 



&L,ti+2 



1 



(|*| 2 ^42,7i+lrfL,7i 



r(l - v4 2 , n+lM, 7i+l) 

+A 1 , n+2 (\r\ 2 - A 2 , n+ iA2,„ + i)d L ,n+2) ■ (A17) 



Eqs. ( [A15| )- flA17p and Eq. (|A8|) above, together with the 
boundary condition that the coefficients should go to 
zero as n goes to plus or minus infinity, determine the 
wave function completely in the region to the left of the 
barrier J13 



Studying Eqs. ( |A16| )-( |A17|) we notice some unneces- 
sary singularities, which we may get rid of by defining 



and 



d L,n = d L ,n/(rrin-lVn+l), 
<7n = 1 - M,nM, n , 



P'n = Vn-lVn+lPn, 
l' n = VnSVn-lln- 



(A18) 



(A19) 



The coefficient b^ is determined from 

^L,77 + 2 = |*| 2 ^42, 71+1-771-1^, 71 

+ J 4i,7i+ 2 (M 2 - A 2tn+1 A 2}n+1 )r) n+3 d' Ltn+2 . (A20) 

Using the primed quantities above there is a new recur- 
sive relation 



4^L,n+2 + P' n d'L,n + 7n rf L,7i-2 ~ Jl$n,t 



(A21) 



which is more suitable to handle from a num erical point 
of view. Solving for d' L we use Eq. ( [AlH ) to get the 
original wave fu nctio n coefficient d^. 

To solve Eq. ( A21 ) we «se a method previously intro- 
duced in the s-wave case.Ea First, we introduce 



, = | d L,J d L,n-2i n > 
rf L,7i/ rf L,7i+2> n < °) 



(A22) 



leading to 



, 1 



a' n x n+2 + P' n + j' n — = 0, n > 

a 'n— + Pn + l'n X n-2 = 0, 71 < 0, 



together with 



(A23) 



(A24) 



Rearranging Eq. (A23) we get a continued fraction ex- 
pansion 



-In 



P'n + a 'n x n+2 



,71 > 0, 



-a r . 



P'n + l'n x r> 



-, n < 0. 



(A25) 



12 



The continued fraction expansion can be evaluated by a 
routine in Ref. ^ giving us a good approximation for 
x ni e , or be truncated at a large n = ni arge . It turns 
out that the truncation method is often very accurate, 
especially if the barrier strength Z is large. Having cal- 
culated x n we get all d' n from the relations 



X2d' L , n > 
x n x n+2 ...x^ 2 d' L . , n < 0. 



(A26) 



There is a similar procedure as above to calculate the 
wave function coefficients of electronlike quasiparticle in- 
cident from the right. One could also get the wave func- 
tions for the leftmovers using the formulas above if the 



-V, 0o 



—4>o and A i 



(2) 



A 



2(1) 



substitutions V 
are made. 

In Eqs. (j|) and (O), the density of states in the super- 
conductors N T {E) is included which will give rise to dis- 
continuities at the gap edges. This will produce numeri- 
cal errors when we perform numerical calculations. For- 
tunately, the transmission amp litude J T appearing in the 
recurrence equation [Eq. ( A15| )] will only enter as a pref- 
actor in the coefficient d, and therefore also in a, b, and 
c. Hence, a factor can be extracted from T T (E, 9, m) 
in Eq. ([lo|) and Eq. ( [l4|) and combined with the density 
of states: 



N T (E)J 2 T (E) 



|f?|-|A r |)(l-^„). (A27) 



The discontinuities at the gap edges are then eliminated. 

Another numerical difficulty is the appearance of a 
number of resonances in the continued fraction expansion 
solution of the recurrence relation. For a particle coming 
in at energy E and angle 9 these resonances appear if 
the trajectory (MAR state) hits zero energy (MGS) or 
the gap edges. When the integration over energy is per- 
formed, it is therefore a good idea to take special care of 
these points in order to faster get an accurate result. 



APPENDIX B: EQUIVALENCE OF CURRENT 
FORMULAS 

In the literature different current formulas have been 
used. In this Appendix we show analytically that they 
are all the same. In addition we have noticed that the 
different formulas numerically give the same results. 

A very useful relation is the completeness relation 

J2 MrR(r') + <(rK(r')] = *(r - r'), (Bl) 



where the sum runs .oyer solutions to the BdG equation 
with positive energy.E2l We apply V r ' and take the imag- 
inary part: 



Y Im {vvVvl) = - ^ Im KVmJ . 



(B2) 



From now on we let the sum over states (here v) run 
over fc-vectors, since we in this paper discuss plane wave 
solutions. 

Next, we have the following symmetry between solu- 
tions to the BdG equation with positive and negative 
energies 



u k {E) 
v k (E) 



-vl k (-E) 



(B3) 



where the left hand side is a solution describing an 
electron-like quasiparticle with energy E, and the right 
hand side is a solution describing a hole-like quasiparti- 
cle with energy — E. Both particles have positive group 
velocities. We can use Eq. (B3) in Eq. (B2) to get the 
same type of completeness also for negative energies. 

The mrrent formula can be written in the excitation 
pictures (positive energies only) as 



3 = — y^y^m-faLVufc* ~v ka Vv* ka } f(E ka ) 
m ■* — ' * — ' 

er k 

eh 



— Y] Y] Im {v k a Vv * ka } 



(B4) 



a k 



where a =f, \ labels spin up and spin down, see 



Fig. |llj(a). Since in this paper a spin independent prob- 
lem is featured, the spin sum just produces a factor of 
two. One of the terms (say a = j) in the spin sum can be 
mapped into negative energies using Eq. (B3) and the re- 
lation f(—E) = 1 — f{E). After this mapping to negative 
energy states the spin sum transforms into a summation 
over branches of negative (/? = —) and positive (/? = +) 
energies (sometimes denoted the semiconductor picture) , 
see Fig. O(b). As a result the current formula after these 
modifications is 



k 



(B5) 



k 



We now study the term in Eq. ( p35| ) without the Fermi- 
factor. We write out the sum over (5: 

^^Imjw^Vw^} 

k 

= H [Im{wfe+Vufc + } -Im{?4_Vufe_}] 

k 

= ^[lm{ Ufc+ V^ + }-Im{^ fe>+ V^ fc)+ }] =0, (B6) 



where we in the first step used Eq. (B2) (in the nega- 
tive energy bransch term) and in the second step used 
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Eq. (B3). We therefore have the following current for- 
mula in the semiconductor picture: 

J ' = ~J2^2 Im i u *kf> yUk f J + V *kl3^V kf 3} f(E kp ). (B7) 

k 

We may now turn the sum over hole-like quasiparticles 
into a sum over electron-like quasiparticles by using the 
symmetry in Eq. (B3) once again. The current formula 
will then be 

J '"^EE Im { U kf3^ u kfJ + vlpVVkp} 



[3 k=k" 



x [2f FD {E kp ) - 1] , 



(B8) 



which is the on e w ritten down in Eq. ([!]). The states 



included in Eq. (|B8|) is shown in Fig. |ll|(c). 

In conc lusion, we can use the completeness relation 
[Eq. (|B2|)1 and the symmetry of the BdG equation [Eq. 



(B3)] to rewrite the formula for the current into different 
forms. Which form to use is a matter of convenience as 
long as one sums over the right set of sc atte ring states. 
For instance in Refs. |l3|, ^J, and [l6] Eq. (B5) was used; 



on the other hand, in this paper and in Rcf. |17j Eq. ([b|) 

has bee: 



1.0 



P 0.5 



0.0 
1.0 



0.5 



-0.5 



-1.0 



. usen 

\ (a)/ 

\ / ° = \ / 


\ / ° = i \ / 


LI 

X/ \ / 

h=\ J \ 
/ (b )\ 





-1.5 
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-0.5 0.5 

k/k„ 
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FIG. 11. Here we show the different sets of states to sum 
over for the different current formulas. In (a) we show the two 
branches with positive energies (labeled by a spin index a) in- 
cluded in Eq. (B4). In (b) we show the negative (/3 = — ) and 
the positive (/3 = +) branches included in Eq. (B7). In (c) we 



show the states included in Eq. (B8): only the electron- like 
particles are included (the solid line), while the hole-like par- 
ticles are excluded (the dashed line). 



APPENDIX C: SOME REMARKS ON THE 
CURRENT 



Letting V — > — V and E — > —E we see from Eq. ( [A3] ) 
that the Andreev reflection amplitude A n —* —A* n . 
Studying Eq. ( A16Q , this gives 



fin *■ fini 
In -> J*, 



(CI) 



implying that Eq. ( A15| ) determining d n changes to 

Q* n d n+2 + l3*d n + j*dn-2 = rJ(E)S nfi . (C2) 

After complex conjugation and multiplication of the fac- 
tor r/r* , we have 



{a n d* n+2 + (3 n d* n + 7„d*_ 2 )— = rJ(E)S n>0 . 



(C3) 



If we compare this recurrence equation to the one in 
Eq. (A15), we see that 



d* n r/r* 



Using Eq. (C4) in Eqs. (AS) and (A17), we get 



* / * 



which implies [from Eq. (|lO[)l that 

T;(E,6,m) -> [TZ(E,0,m)]\ 

A minus sign is extracted when we let E 
Eq. (||) since 

tanh(-.E/2£; B T) -> -tanh(-E/2k B T). 



(C4) 



(C5) 
(C6) 
(C7) 



(C8) 
-E in 

(C9) 



Eq. 



Taking Eq. ( ]C8[ ) and Eq. ( |C9[ ) into account we get from 



.4 



x,0 



-A 



x,0- 



Therefore, one finally finds 

i*A-v) = -£.,o(V). 



(C10) 



(Cll) 



The same kind of reasoning as above for the pe rpen - 
dicu lar ca se is used to put the transformations Eqs. ( |C4| )- 
dC7|) and (C9) into the corresponding equations for j y Q . 
Therefore, we also have 



JyA-V) = -Jy, (V). 



(C12) 



1. Derivation of the symmetry /(— V) ~ —I(V) 



In this subsection we show the relation I(—V) = 
—I(V) for both perpendicular and parallel current den- 
sity. 
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2. Disappearance of ac components 

Here we set out to show analytically how the ac 
components disappear for some orientations of the su- 
perconducting order parameters. E3 The argument pre- 
sented here is not ^-dependent: it is true for any barrier 
strength. 

The keypoint is to compare current contributions as- 
sociated with injection angles ±9 and see how they add 
up. We can separate the integrand which we integrate 
over the angles 9e[— n/2, tt/2] into two functions: first a 
trigonometric function and second a complicated func- 
tion depending on the gaps. The trigonometric function 
is cos 9 (even in 9) for the x current and sin 9 (odd in 
9) for the y current. We will show that the second func- 
tion is odd or even in 9 depending on the orientation of 
the superconductors. The combination of the two func- 
tions will therefore be odd in some cases and therefore 
integrate to zero. 

Case 1: do\do junction. In this case the gaps are sym- 
metric in 9 on both the left and the right side, implying 
that the second function is even in 9. Taking the trigono- 
metric functions into account, the x current therefore 
contains all components m while the y current is zero for 
all m. 

Case 2: junction. This case is shown in Fig. |^. 

The gap on the right side is oriented so that the negative 
lobe appears for negative angles: A2(— 9) = —A2(9) = 
A2(9) ex-p(iir). This extra phase ir for the negative angles 
can be given to c/)q. As noted in Appendix |A|, the phase 
difference over the junction (j)Q is factored out in the be- 
ginning of the calculation, but it reappears in the phase 
fl m in the final expression for the current, see Eq. (11). 
We therefore have 



n m (-6) = Q, m {6) +TO7T. 



(C13) 



In addition, we have (taking also the trigonometric func- 
tion into account) 



A x , m {~ 9) = A xm (9) 
B x , m (—9) = B XyTn (6) 

Ay im { 9) = ^^y,m{9) 
, B y , m ( — 0) = -By tm {9) 



Oi x ,m{-9) = Ci x . m {9) 

Cy : m( — 9) = —Cy i7n (9) 

(C14) 



Using Eq. (C13) and Eq. (C14), we can rewrite the x 
current 



- = VU«>o)(i 

and the y current 

- = £<W0>oxi 



)e 



:(a„ 



») 



Jim 7T\ Jl 



(C15) 



(C16) 



where the functions C x / ViTn (9 > 0) contains an integra- 
tion over positive angles only. It is clear from Eq. (|C15|) 



[Eq. ( |C16 )] that the odd (even) components of the x (y) 
current will be zero for this orientation. Note that the 
m = component of the y current is also zero. 

Case 3: d^/^d^/^ junction. Here we have an extra 
phase 7r on both the left and the right side for negative 
angles. This means that we have no effective extra phase 
that we can factor out together with (f>Q. Instead we see 
in Eq. (A3) that the Andreev reflection amplitudes are 
A/2(-9) = -A 1/2 {9) = A 1/2 (6>)exp(i7r). Studying the 
equations in Appendix [a], only combinations of two A 
multiplying each other appear. Thus the extra phases 
7r cancel. Taking also the trigonometric function into 
account, the integrand will therefore be even (odd) for 
all m components of the x (y) current in the same way 
as for the do\do junction. 

For orientations in between ai/ 2 = and ai/2 = 7r/4, 
the cancellation effects will not be perfect and all com- 
ponents m are non-zero for both the x and y currents. 



APPENDIX D: ANALYTIC EXPRESSIONS FOR 
LARGE AND SMALL VOLTAGES 

1. Large voltage in the ballistic limit 

In the large voltage limit (eV » A ) we only need to 
consider one Andreev reflection. This means that only 
coefficients a n , b n , c„, and d n with indices n — and 
n = ±2 have to be taken into account. 

In the ballistic limit (Z = 0), the surviving coefficients 
from an electron-like quasiparticle injected from the left 
superconductor at energy E and angle 9 will be a^* (the 
injected electron) and b^* (the hole created from Andreev 
reflection at the right superconductor). As outlined in 
Ref. ^ we can then perform the integration over energy 
analytically in the limit eV — > oo and get the following 
expressions for the currents 



(TO 



3v,o 



(CTO/Ly 



eV 


2 


Ao 


+ 3 


eV 






+ Ie 


Ao 






• tt/2 


u 


-tt/2 



tt/2 



d9 cos 9 



-tt/2 



d9 sin ( 



|Ai(0)| , |A 2 (0)| 



A 



Ao 



IM0)| , |A 2 (0)| 



(Dl) 



valid for large voltages. Note that both the excess current 
and the surface current is independent of voltage when 
eV » A . 



2. Small voltage in the ballistic limit 

In the ballistic case there is a nonzero current contri- 
bution in the small-voltage limit for both the perpen- 
dicular (x) and parallel (y) direction. For the ballistic 
one-channel SNS junction with different s-wave gaps one 
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has previously derived that the perpendicular current is 
proportional to the smallest of the gaps.CJ Generalizing 
this result to our case we have for the current in the 
x-direction 

^cos.min^,^). (D2) 

-tt/2 An An 



k 

(TO 



In the same way we find for the parallel current density 
in the normal region that 



Jy 



(<Ja/L y ) 



r/2 



d6 sin #min( 



ia 2 (0)i 



-tt/2 



An 



(D3) 
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